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Recent atomic resolution scanning tunneling microscope (STM) measurements discovered remark- 
able electronic inhomogeneity, i.e. nano-scale spatial variations of the local density of states (LDOS) 
and the superconducting energy gap, in the high-T c superconductor Bi2Sr2CaCu2 08+a;. Based on 
the experimental findings, we conjectured that the inhomogeneity arises from variations in local oxy- 
gen doping level and may be generic of doped Mott insulators. In this paper, we provide theoretical 
support for this picture. We study a doped Mott insulator within a generalized t-J model, where dop- 
ing is accompanied by ionic Coulomb potentials centered in the BiO plane located a distance d s away 
from the Cu02 plane. We solve, at the mean-field level, a set of spatially unrestricted Bogoliubov-de 
Gennes equations self-consistently to obtain the distributions of the hole concentration, the valence 
bond and the pairing order parameters for different nominal/average doping concentrations. We 
calculate the LDOS spectrum, the integrated LDOS, and the local superconducting gap as those 
measured by STM, make detailed comparisons to experiments, and find remarkable agreement with 
the experimental data. We emphasize the unconventional screening of the ionic potential in a doped 
Mott insulator and show that nonlinear screening dominates on nano-meter scales, comparable to 
the short coherence length of the superconductor, which is the origin of the electronic inhomogeneity. 
It leads to strong inhomogeneous redistribution of the local hole density and promotes the notion of 
local doping concentration (LDC). We find that the inhomogeneity structure manifests itself at all 
energy scales in the STM tunneling differential conductance, and elucidate the similarity and the 
differences between the data obtained in the constant tunneling current mode and the same data 
normalized to reflect constant tip-to-sample distance. We also discuss the underdoped case where 
nonlinear screening of the ionic potential turns the spatial electronic structure into a percolative 
mixture of patches with smaller pairing gaps embedded in a background with larger gaps to single 
particle excitations. 



PACS numbers: 74.20.Mn, 74.80.-g, 68.37.Ef, 74.25. Jb 



I. INTRODUCTION 

The parent compounds of the high-T c superconduc- 
tors are stoichiometric Mott insulators pi. In principle, 
carriers can be introduced into the insulating state ei- 
ther by field effect or by chemical doping, but to date 
high-r c superconductivity in the cuprates has only been 
found to arise when they are properly doped away from 
stoichiometry. In the case of YBa2Cu307_ a; (YBCO) 
and Bi2Sr2CaCu208+a: (BSCCO), the dopants are the 
excess oxygen atoms. The latter inject holes into the 
CuC"2 planes, leaving behind the negatively charged oxy- 
gen ions. While for YBCO, the dopant ions go into the 
copper oxide chains and can be arranged more period- 
ically, M in the case of BSCCO, the dopants, with a 
concentration of x, disorder themselves in the semicon- 
ducting BiO layer which is about 5A away from the Cu02 
plane where superconductivity is believed to originate. 

In addition to doping holes into the Cu02 planes, 
the dopant oxygen atoms inevitably introduce long-range 
random ionic potentials that scatter the carriers in the 
copper-oxide plane in analogy to the situation encoun- 
tered in two-dimensional (2D) electron systems in mod- 
ulation doped semiconductors. This has the potential to 



cause the electronic structure in the copper-oxide plane 
to be inhomogeneous. Early neutron scattering |jj, tun- 
neling H and STM || data showed features that could 
be accounted for by inhomogeneity, but they were usually 
attributed to sample quality. The common view, which is 
largely based on the physics of ordinary metals, has been 
that the ionic potential in a high-quality crystal would 
be screened by phonons and by the carriers in the plane 
such that charges distribute uniformly and are accommo- 
dated by a homogeneous electronic structure; albeit that 
the latter has the tendency towards microscopic phase 
separation H due to strong electronic correlations. This 
viewpoint appears to be supported, to a certain extent, 
by transport M and photoemission S measurements on 
BSCCO which observe quasiparticles with relatively well- 
defined momentum, a mean-free path of about 100A in 
the nodal direction, and a transport mean-free path of 
about 400 A. 

However, this conventional view has been seriously 
challenged by recent low-temperature STM measure- 
ments on BSCCO. @-0 In Ref. §, Pan et al. (here 
after referred to as we) observed spatially cross-correlated 
variations in the LDOS and the superconducting energy 
gap on a remarkably short length scale of about 14A. 
We identified this as the gap amplitude coherence length 



or the pair size. It directly confirms the local pairing 
nature in the high-T c superconductor. The magnitude 
of the superconducting gap has a Gaussian distribution 
with a mean value of 42 meV, close to the gap value ob- 
tained from planar tunneling, |l2| and a full-width-half- 
maximum of about 20 meV. This variance coincides with 
the intrinsic width (~ 20 meV) of the coherence peak 
measured by angle resolved photoemission spectroscopic 
(ARPES) near the antinodes where the d-wave gap is at 
its maximum. |13|] It offers an explanation of the latter 
as arising, even in the absence of bilayer splitting, from 
averaging over a distribution of local gap values under 
the macroscopic ARPES light spot. 

Careful analysis of the STM data, combined with those 
of the ARPES, led us to conjecture in Ref. || that the 
observed inhomogeneous electronic structure arises from 
the ionic potential associated with the off-stoichiometry 
oxygen dopants disordered in the BiO layer. We argued 
that the screening of the ionic potential in the doped 
Mott insulator is unconventional and incomplete on the 
nano-meter scale of the pair size, causing the local doping 
(hole) concentration to exhibit spatial variations. Such a 
doping profile results in a highly cross-correlated pairing 
gap profile if the pairing interaction is sufficiently short- 
ranged. 

It should be stressed that the physics involved here 
is highly unconventional. In an ordinary metal such as 
copper, the Thomas-Fermi screening length in the lin- 
ear screening theory is shorter than lA, beyond which 
the ionic potential is perfectly screened and the charge 
density homogeneous. Such a short screening length is 
a result of the large thermodynamic density of states or 
the compressibility (dn/dfi) in ordinary metals. However, 
adding holes to a doped Mott insulator causes significant 
shift in the chemical potential due to the strong local 
Coulomb repulsion. This leads to a small electronic com- 
pressibility and thus a larger screening length even within 
the linear screening theory. We believe that the observed 
inhomogeneous superconducting state is a result of the 
strong local correlation in a doped Mott insulator. 

It is important to provide theoretical tests for these 
ideas and to compare theoretical results with experi- 
ments. In this paper, we do so by studying the d-wave 
superconducting state of a doped Mott insulator in a mi- 
croscopic theory of a generalized t-J model, where dop- 
ing is accompanied by ionic Coulomb potentials centered 
a distance d s away from the 2D plane. Since an exact 
solution of this model for large systems is not possible 
even numerically at present, we provide a self-consistent 
mean- field solution that captures the essential physics. 
The hope is that upon comparing the results to exper- 
iment, we can develop useful insights to guide our un- 
derstanding in the interim. Specifically, we solve a set 
of spatially unrestricted Bogoliubov-de Gennes equations 
self-consistently at every site on square lattices of up to 
32 x 32 sites to obtain the spatial distributions of the 
hole concentration, the resonating valence bond and the 
pairing order parameters for different nominal/average 



doping concentrations. We find that the inhomogeneity 
originates from nonlinear screening of the ionic potential. 
We show that a single negatively charged test ion inserted 
at a setback distance d s above the 2D plane induces a 
nonlinear screening cloud wherein the doping concentra- 
tion is significantly larger than its averaged value. For 
a given strength of the Coulomb interaction, the size of 
the nonlinear screening cloud is controlled by the distance 
d s and is independent of the average doping concentra- 
tion. With a finite density of dopant ions (same as the 
averaged doped hole density), an inhomogeneous elec- 
tronic structure on the scale of d s emerges that shares an 
analogy to the 2D electron system in modulation doped 
semiconductors where nonlinear screening leads to an in- 
homogeneous mixture of metallic and dielectric regions 
at low electron density. |14| The immediate consequence 
of the nonlinear screening is the strong spatially inho- 
mogeneous redistribution of the charge carrier density 
brought about by the bare ionic potential. The concept 
of a local doping concentration (LDC) || therefore natu- 
rally arises. The d-wave pairing amplitude turns out to 
follow locally the spatially varying hole density, resulting 
in an inhomogeneous superconducting state. In the res- 
onating valence bond picture M, this corresponds to an 
inhomogeneous state of spinon pairing and local holon 
condensation. We calculate the local tunneling density 
of states spectrum, the integrated LDOS, and the local 
superconducting gap as those measured by STM. Sur- 
prisingly, the local spectral properties can be well de- 
scribed by the Mott-Hubbard picture once the notion of 
a LDC is established. We perform a detailed comparison 
of our results with the experimental data and find re- 
markable agreement in support of the conjectures made 
in Ref. [B| and the theoretical picture of spinon pairing 
with local holon condensation for the inhomogeneous su- 
perconducting state. Extending this picture to the un- 
derdoped regime, where the averaged inter-hole separa- 
tion becomes larger than the setback distance d Sl nat- 
urally leads to percolative structures of superconducting 
patches immersed in a background with a large tunneling 
gap to single particle excitations. 

Currently, there exists a "discrepancy" in the pub- 
lished STM data by different groups p-|ll|| with regard to 
whether the LDOS shows spatial inhomogeneity at low 
energies. We point out that the difference in the pre- 
sented experimental data depends on whether they have 
been normalized to remove certain matrix element ef- 
fects. In Ref. H , we have taken into account, by correctly 
normalizing the data, the effect of the tunneling matrix 
element along the direction perpendicular to the surface. 
The normalized data reflect the tunneling spectra when 
the tip to sample distance is kept effectively at a con- 
stant. We believe this is physical. The low energy elec- 
tronic spectrum as measured by the local tunneling den- 
sity of states is indeed inhomogeneous. H In this paper, 
we show that the same inhomogeneity structure mani- 
fests itself at all energy scales in the calculated tunneling 
spectra. We find that the calculated zero-bias tunneling 



conductance, which is dominated by the nodal quasipar- 
ticles, shows spatial variations that are correlated with 
that of the LDC. We elucidate in detail the mapping, the 
differences and similarities between the tunneling differ- 
ential conductance obtained in the constant tunneling 
current mode, which docs not directly represent the elec- 
tron LDOS, and the same data normalized to reflect a 
constant tip-to-sample distance that faithfully represent 
the electronic tunneling density of states. We also pro- 
vide preliminary evidence that the high superconducting 
transition temperature T c can be protected by the short 
superconducting coherence length and coexist with this 
type of electronic inhomogeneity. 

The rest of the paper is organized as follows. In sec- 
tion If, we describe our generalized t-J model and set up 
the self-consistent equations for the spatially unrestricted 
Hartree-Fock-Bogoliubov mean-field solutions. In section 
III, we study the screening properties of the ionic poten- 
tial in detail. A single test ion is inserted which imposes a 
coulombic potential on the otherwise uniform d-wave su- 
perconducting state. The charge redistribution and the 
response of the pairing order parameter to screening of 
the ionic potential is analyzed. In section IV, we discuss 
the solution of the inhomogeneous d-wave superconduct- 
ing state for a finite density of dopant ions and doped 
holes. We emphasize the concept of LDC and analyze 
the distribution of the d-wave order parameter and its 
correlation with the local hole concentration in subsec- 
tion IV. A. The numerical results of the local tunneling 
spectrum are presented in subsection IV. B and compared 
to experimental data. We analyze the statistical proper- 
ties of the LDC, the integrated LDOS and the pairing 
gap distribution and the spatial correlation and cross- 
correlation among these physical observables. In section 
V, we focus on the inhomogeneity at low energy scales. 
We study the spatial distribution of the zero-bias tunnel- 
ing conductance, its correlation with the LDC, and make 
predictions for experimental tests. In section VI, we dis- 
cuss the mapping between constant-current differential 
tunneling conductance and the constant tip-to-sample 
distance LDOS spectra. The STM topography is calcu- 
lated for our system together with the constant-current 
differential tunneling conductance and compared to ex- 
periments. A summary of the results and their implica- 
tions are given in section VII together with discussions 
of several open issues. 



II. GENERALIZED t J MODEL AND THE 
SELF-CONSISTENCY EQUATIONS 

We begin with the generalized t-J model Hamiltonian 
that includes the long-range Coulomb interaction and the 
ionic potentials introduced by the process of carrier dop- 



ing, 



H = H t -j + H t 



Coul 



Hu 



(1) 



Here H t -j is the usual 2D t-J model on a square lattice, 



(id) 



(id) 



where c ia is the electron creation operator and Si is the 
spin operator Si = {\)c\ a a a pCij3. The sums over (i,j) 
are among nearest neighbors and sums over repeated spin 
indices are implied. The most important Mott-Hubbard 
physics, i.e. the strong on-site Coulomb repulsion, is in- 
cluded in the additional constraint of no double occu- 
pancy at each site, rii — c\ a Ci a < 1. Hq ou \ is the long- 
range Coulomb repulsion, 



H, 



Coul 



Y,V t n t , Vi = V c 



E 

3^ 



(3) 



where n denotes the average density. The Coulomb in- 
teraction strength V c is given by V c — e 2 /4irea, with 
e ~ 8 the dielectric constant (jj| and a ~ 3.8 A the lat- 
tice constant corresponding to the Cu-Cu atomic spacing. 
The resulting V c ~ t. As emphasized in Ref. |l(|, the 
inclusion of the long-range Coulomb interaction is nec- 
essary to prevent the mean-field ground state of H t -j 
from macroscopic phase separation in the interesting pa- 
rameter regime. |Ij] The off-plane ionic potential of the 
dopants is described by -ffi on , 



Ni, 



H u 



y^ Uim, Uj = y^ 



v io 



=i v>r 



di 



(4) 



Here d s is the distance between the C11O2 plane and the 
BiO layer where the negatively charged ions reside ran- 
domly at ri on , and N- lon is the number of independent 
ions in the BiO layer. To model the situation in BSCCO, 
where each dopant oxygen gives one hole to each of the 
planes in a bi-layer, we use A^ ion = iVhole = x x N s where 
AWc is the number of doped holes, x is the average dop- 
ing (hole) concentration on a lattice of N s sites. We use 
J as the unit of energy and set d s = 1.5a, V{ on — V c = 5 J, 
and t = 3 J in most of our numerical calculations. We ver- 
ified that varying these parameter in a reasonable range 
does not qualitatively change our results. 

It is convenient to describe the projected Hilbert space 
in terms of a spin-carrying fermion, the spinon fj a , cre- 
ating the singly occupied site with spin-er and a spinlcss 
boson, the holon bi, keeping track of the empty site p8). 
The electron creation operator becomes c ia = f ia bi and 
the occupancy constraint translates in this slave-boson 
formulation into f icr fia + b\bi = 1. In the mean-field 
theory, the antiferromagnetic spin-exchange term is de- 
coupled according to |fj| 

3 r 



Oj • 3, 



M j fUj<r + ^ 



+ A* J (f li f^-f i] f n )+h.c. 
+ o(lxd 2 



|Ay| 2 ), 



where Ay and Xij are the spinon pairing and bond order 
parameters respectively, 



(/a/jT ~ Mfji): Xij - {fiafjo), 



(5) 



defined for each nearest neighbor bond. The inhomoge- 
neous superconducting phase is reached through the local 
condensation of bosons at low temperatures, 



(M) = (bi) 



(6) 



It is important to emphasize that in the presence of the 
translation symmetry breaking ionic potential in Eq. (H), 
(A^, Xijibi) become spatially dependent and must be de- 
termined at every site self-consistently. Note that the lo- 
cal doped hole density or the LDC is then directly related 
to the local boson condensate density, 



= 6 



(7) 



and the constraint at the mean-field level becomes nj — 

flafi<? = 1 — x ii which is enforced on average at every 
site by locally shifting the fcrmion chemical potential /i/ 
to Xi + fij. Throughout this paper, we use x to denote 
the average (doped) hole density or the nominal doping 
concentration, 



1 N s 

— Y 



(8) 



The decoupled Hamiltonian can be written down using 
the Bogoliubov-Nambu formalism over the Hilbert space 
of paired spinons, 



# = £(#,/< 



U) 



K, 
FT. 



F ■ 
-K* 



/it 

rt 



f] 



(9) 



where the sums on i and j run over all lattice sites and 






(10) 



Ku = -{tb\ + -J Xij ) J2 hi+v + [V sc (i) - HflSij, (11) 



(12) 



with 7/ = ±x, ±y. In Eq. (|l l[) , V sc (i) given by 
V ac (i) = U t + A, + V c V J '~ g 



jV« 



is the nonlinearly screened local potential seen by the 
spinons implied by the self-consistency conditions. Note 
that in deriving Eq. (U2), only the Hartree potential in 
the long-range Coulomb interaction is retained, whereas 
the exchange potential is neglected. The effects of the 
latter on the superconducting state will be studied else- 
where. 



The Hamiltonian in Eq. (||) can be diagonalized in 
real space by solving the corresponding Bogoliubov-de 
Gennes equations to obtain the eigenstates 7^ and j n 
with energy E n ,n — 1, ...,2N S . The spinon operator 
can be expanded in this basis according to 



/J(*) = E u »w^ e_iBB 



L/Ti 



filV) = £ w ™(«)7^ 



-iE„t/h 



(13) 
(14) 



where (u n (i) , v n (i)) is the wave-function at site i. The 
order parameters and the local hole density can be ex- 
pressed in terms of the wave-functions as follows 

Ay=£[<(*K0')[l-/(2k)] 

-u n {i)v* n (j)f{E n % (15) 

x« = ZK(*xo')[i-/(^)] 

■II 

+ u n (t)u*Jj)f(E n )], (16) 

l-x i = '£[\v n (i)\ 2 [l-f(E n )} 

II 

+ \u n (i)\ 2 f(E n )], (17) 

where f(E n ) is the usual Fermi distribution function. 

We solve the self-consistency Eqs. (pfjl7|) through nu- 
merical iterations. Typically, we start with a random 
set of (Aij,Xij, K,%i, Hf), insert them into Eq. (0) and 
diagonalize the resulting matrix to obtain the wave- 
functions (u n (i) , v n (i)) and the eigenvalues E n . Then 
we update the set of (Aij,Xiji^i,Xi,iJ,f) according to 
Eqs. (|l5||lcUl7| ) and insert them back into the Hamilto- 
nian (M). The procedure is iterated until convergence is 
reached. In general we allow Xij ancl A,j to be complex 
at the beginning of the iteration, but both of these or- 
der parameters converge to real values at the end of the 
iteration in the average doping range (0.06 < x < 0.32) 
studied. 

In order to compare to STM data, we calculate the 
local tunneling density of states, 

N i (u)=lm[dte iut G% t (t), (18) 

where GJJ* (t) is the retarded local Green's function of the 
electrons, 



GT{t) = -e(t){{cl(Q),c l .{t)}). 



(19) 



In the mean-field theory, since the bosons locally con- 
dense, we have 



Gr(t)=-W)<{&(0),/.v(t)}> 



(20) 



Notice that the local electron Green's function and there- 
fore the local tunneling density of states is equal to the 



local holon condensate density times the local spinon 
Green's function, PCJ From Eq. (M), it is clear then that 
the LDOS is overall proportional to the LDC Xi. Using 
Eqs. (|l|,0), we obtain, 

N^lo) = x^ [\u n {i)\ 2 6(tu - E n ) + \v n (i)\ 2 S(uj + En)] . 

H 

(21) 

The role played by the LDC is two-fold: it enters as an 
overall prefactor through essentially the wave-function 
renormalization of the electrons as well as through the 
LDOS of spinons from the response of Ay and Xij to a 
spatially varying Xi . This is a property of the tunneling 
density of states which is different from the thermody- 
namic density of states. pifl 

Our choice of the t-J model for quantitative calcula- 
tions is a natural one, because it is the simplest model 
that captures the physics of a doped (antiferromagnetic) 
Mott-insulator. In this strong coupling approach, the 
carrier concentration is proportional to the doping con- 
centration (rather than the electron density) , which is the 
most fundamental property of a doped Mott-insulator. 
More specifically, as can be clearly seen from Eq. (J2(j) , the 
coherent weight of the quasiparticle is small at small dop- 
ing and scales with the LDC. This is in excellent agree- 
ment with the recent ARPES measurements [|13],£2| on 
BSCCO, where the weight of the emergent quasiparticle 
peak below T c near the anti-node was found to be propor- 
tional to the doping concentration at low temperatures. 



III. A SINGLE OFF-PLANE ION AND 
NONLINEAR SCREENING 

To understand the screening properties of the ionic po- 
tential better, we start with the case where a single off- 
plane ion (i.e. N- lon = 1) is placed above the center of a 
25 x 25 lattice with a setback distance d s — 1.5a. The 
ion imposes a Coulomb potential on the otherwise uni- 
form superconducting state at a doping level x. This 
situation is similar to that of a single nonmagnetic im- 
purity (say Zn) extensively studied using the t-J model 
recently, [|23|-|25| except that the impurity is located out 
of the plane and the impurity potential is long-ranged. 
The scattering is presumably not in the unitary limit. 
In addition, fully self-consistent solution in (pi, Ay, Xij) 
has not been studied before. 

Fig. 1 shows the bare ionic potential Ui from Eq. (0) 
and the screened potential V sc (i) from Eq. jl^ ) coming 
out of the self-consistent solution. The average doping 
concentration in the plane is x = 0.12. It is clear that 
the long-range impurity potential in Fig. la is perfectly 
screened at distances much larger than d s but poorly 
screened at short distances resulting in a short-range im- 
purity potential in Fig. lb. 




FIG. 1. The Coulomb potential from a single negatively 
charged ion located at a distance d s = 1.5a above the 
center of an otherwise uniform 2D d-wave superconductor. 
The lattice size is 25 x 25 and the average doping x = 0.12. 
(a) The long-range bare ionic potential from Eq. (0). (b) 
The screened potential from Eq. ([L2|). 
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FIG. 2. Spatial variations of the LDC induced by the 
ionic potential in Fig. 1. (a) The 2D plot of the local 
hole density (xi) redistribution as a result of nonlinear 
screening, (b) The hole density and the corresponding 
d-wave pairing order parameter (OP) along a horizontal 
line-cut passing through the center of the 25 x 25 lattice at 
x = 0.12. (c) The azimuthally averaged auto-correlation 
function of the spatial variations in local hole concentra- 
tion 8xi defined in Eq. (E2j). The solid line is a Gaussian 
fit to the data, which gives a decay length £ CT = 1.86a. 

The nature of the screening is revealed when we look 
at the redistribution of the carrier density as a result 
of screening. The self-consistently determined local hole 
concentration Xi is plotted in a 2D hole density map in 
Fig. 2a, which shows that in the vicinity of the dopant 
ion, the hole density is strongly inhomogeneous, namely 
the changes in the hole concentration Sxi , brought about 
by the dopant ionic potential, are comparable to the av- 
eraged concentration x itself. Fig. 2b shows the hole 
density distribution along a horizontal line-cut passing 
through the center of the lattice. This clearly demon- 
strates that the screening of the ionic impurity potential 



in the doped Mott insulator is highly unconventional and 
is entirely dominated by nonlinear screening J14J on the 
scale of the ionic setback distance d s . The nonlinearality 
in screening is further verified by its nonlinear response 
to changes in the strength of the Coulomb potential V\ on . 
It is instructive to compare to conventional screening 
in the linear response theory described by the Thomas- 
Fermi screening length, which is applicable to ordinary 
metals and superconductors. For an electron gas in 
2D, the Thomas-Fermi screening length is independent 
of electron density and is given by q^ F = a* B /2, where 
a* B = h 2 e/m*e 2 is the effective Bohr radius. |g6[ Taking 
e ~ 8 and a thermodynamic effective mass m * /m ~ 3, 
which is appropriate for the cuprates, gives a* B ~ 1.4A. 
The resulting Thomas- Fermi screening length ~ 0.7 A 
which is much less than a lattice spacing a ~ 3.8A Thus 
linear screening theory would have predicted that the 
ionic Coulomb potential is screened on the length scale 
of a lattice spacing and the charge redistribution induced 
by the bare potential is small. In our case, Figs. 2a and 
2b clearly show that linear screening is only valid when 
the distance to the center of the ionic potential is much 
greater than the distance d s , where the hole density is 
essentially uniform and its value close to the average con- 
centration x — 0.12. 

Thus, the profile of the screened impurity potential 
and that of the hole density distribution and, as we shall 
see later, that of the d-wave pairing OP (see Fig. 2b), 
are controlled by the crossover from nonlinear to linear 
screening as one moves away from the center of the ionic 
potential. To better quantify this crossover behavior, we 
study the auto-correlation function of the local hole den- 
sity variations 5xi = Xi — x, 



C x (r 3 ) 



— ^2{5xi5x i+ j) 



(22) 



In Fig. 2c, we show the azimuthally averaged C x , which 
decays very fast with the distance, indicative of strong 
short-range correlation in the LDC. The short-distance 
behavior can be fitted very well by a Gaussian with a de- 
cay length £ a = 1.86a. For convenience, let us define the 
correlation length by the decay length of the Gaussian 
function. We arrive at a hole-density correlation length 
£ x = 1.86a, which is quite close to the setback distance 
d s = 1.5a, suggesting the latter as the length scale over 
which the crossover between linear and nonlinear screen- 
ing takes place. 
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FIG. 3. The average doping and setback distance de- 
pendence of screening. (a) The spatial variations in 
hole concentration (Sxi) along the central line-cut at four 
different average doping concentrations, showing screen- 
ing is insensitive to average doping x. (b) The decay 
lengths of the hole-density (diamonds) and the d-wave 
pairing OP (triangle) auto-correlations, and that of the 
cross-correlation (circle) between the two as a function of 
average doping x. (c) Same as in (a), but for a larger 
setback distance d s — 5a. 

We next study the average doping dependence of the 
screening in the nonlinear regime. In the linear screen- 
ing theory, the fact that the 2D Thomas-Fermi screening 
length is independent of carrier density is a consequence 
of the properties of the density of states and the polariz- 
ability of the 2D electron system. |26| This suggests that 
the insensitivity of screening to carrier concentration in 
2D should be true even for nonlinear screening. Analyti- 
cal solutions of the 2D nonlinear screening problem only 
exists for the special case of an antidot in a semicon- 
ductor 2D electron gas in the continuum, gj To study 
this question in the d-wave superconducting state of a 
doped Mott insulator in the presence of a periodic lat- 
tice potential, we plot, in Fig. 3a, the variation of the 
local hole density from its average value (5xi) as a func- 



tion of distance along the center horizontal cut for several 
different average hole concentrations. The fact that the 
data points almost completely collapse onto a single curve 
shows that the screening is indeed insensitive to the car- 
rier density in both the linear and the nonlinear screening 
regime. This is verified quantitatively by the very weak 
doping dependence of the decay lengths extracted from 
the Gaussian fits to the correlation functions at the cor- 
responding doping levels, shown in Fig. 3c. Thus, for a 
given strength of the ionic Coulomb potential, the size 
of the nonlinear screening cloud around the center of the 
potential is controlled entirely by the setback distance d s . 
In Fig. 3c, we show the spatial variation of the LDC 8xi 
along the center horizontal cut for the case of a larger 
setback distance d s = 5a at x = 0.06 and 0.12. It is clear 
that as the screened potential becomes smoother, the size 
of the screening cloud becomes larger, characterized by a 
decay length on the order of d s . 
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FIG. 4. The 2D maps showing the spatial variations 
in the dominant d-wave (a), and the accompanying ex- 
tended s-wave (b) pairing OPs in the vicinity of the ionic 
potential. 

In Fig. 4, we show the spatial variation of the self- 
consistently determined superconducting order parame- 
ter near the center of the ionic potential for d s ~ 1.5a. 
We find that the nonlinear screening region where the 
hole concentration exhibits strong inhomogeneity is ac- 
companied by spatial variation of the order parameter 
Aij. Moreover. 



when A^, defined in Eq. (||), is decom- 



posed into a d a .2_ J; 2 and an extended s-wave component 
at a site i according to, 

A-d{i) = -7\^is+x + ^is-x ~~ A iii+ y — Aii_y), (23) 
A s (i) = -(A M+£ + A M _ £ + A t , l+y + A M _ 9 ), (24) 

we find that the spatially varying d-wave order parame- 
ter Ad(i) shown in Fig. 4a is complimented by a much 
smaller nonuniform A s (i) shown in Fig. 4b. The magni- 
tude of Ad along the center line cut in Fig. 4a has been 
shown in Fig. 2b. 

This situation is quite reminiscent of a single in-plane 
nonmagnetic impurity, e.g. an intentionally doped Zn 
atom, that replaces a copper atom causing strong scatter- 



ing in the unitary limit. |23|- |25|J28| , |29| The inhomogene- 
ity induced A s shows an interesting four-fold (d-wave) 
symmetry, i.e. it is vanishingly small along the nodal di- 
rection of Arf(i) and changes sign upon a 90° rotation. 
These properties can be qualitatively understood from 
the Ginzburg-Landau theory |23|,|30| that permits by sym- 
metry a mixed gradient term proportional to dAddA s . 
Note that the order parameters in this d + s state are 
both real and that spontaneous time-reversal symmetry 
breaking, such as in a d + is state, does not occur. It is 
interesting to point out that the reduction is greater in 
the anti-nodal direction where the d-wave gap is at its 
maximum. 
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FIG. 5. The azimuthally averaged auto-correlation 
function for the d-wave OP variation (a) and the abso- 
lute value of the cross-correlation function between the 
hole density variation and the d-wave OP variation (b). 
The solid lines are Gaussian fits to the data. 

The auto-correlation function of the spatial variation 
of the d-wave OP SAd(i) = Ad(i) — Ad and the cross- 
correlation function between 5Ad(i) and the local doping 
variation Sxi are defined by 



C x -A d { V J 



= j^Y,( SA ^ SA ^ +j))> 



i 



Z^ 



(5xi6A d (i+j)). 



(25) 
(26) 



We plot the azimuthally averaged correlation functions in 
Fig. 5a and 5b. The A^-correlation length ^A d — 2.53a, 
which is close to the hole density correlation length. The 
fact that the cross-correlation length, £ X -A d ~ 2.26a, is 
close to £ x and £,A d shows that the spatial variations in 
Ad and the LDC are strongly (anti-)correlated. This is an 
important character of the d-wave superconducting state 
emerged from the short-range resonating valence bond 
state, |y|yj where the local spinon pairing produced by 
the nearest-neighbor antiferromagnetic exchange inter- 
action is directly affected by the local holon condensate 
density. 

To summarize this section, an out-of-plane negatively 
charged ion induces a nonlinear screening cloud in its 
vicinity, wherein the doping concentration is significantly 
higher and the d-wave pairing parameter significantly 
smaller than their averaged values in the bulk of the 



2D superconductor. The size of the nonlinear screening 
cloud is essentially independent of the averaged doping 
concentration. It is controlled by the setback distance d s 
of the dopant ions to the 2D plane. This local picture 
for the response of the d-wave superconducting state to 
a single ionic test charge will be important for under- 
standing the real situations of many ions, particularly in 
the case of dilute dopants in the underdoped regime. In 
the rest of the paper, we use d s = 1.5a. This comes out 
naturally from the lattice constants of BSCCO which is 
about 3.8A in the plane and 5A in the c-direction. We 
will come back to this issue in the last section and discuss 
other physical processes, such as dielectric screening from 
interband transitions and electron-phonon coupling, that 
could complicate the choice of d s , but would not change 
the essential physics described in this paper. 
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FIG. 6. 2D density plot of the bare (a) and the screened 
(b) ionic potential on a 32 x 32 lattice at average doping 
x — 0.12. The data has been smoothed to reduce the 
effects of the lattice discreteness. 



IV. FINITE DENSITY OF IONS — 

INHOMOGENEOUS SUPERCONDUCTING 

STATE OF A DOPED MOTT INSULATOR 

Having studied how an external out-of-plane charge 
is screened by the d-wave superconductor, we now turn 
to the situation in BSCCO where iVi on number of neg- 
atively charged ions resides in the BiO layer donating 
-Whole = -Wion number of doped holes to each of the CuO 
bilayers. Without the knowledge of the detailed chemical 
frustration encountered by the dopant ions, we assume 
that they distribute randomly in the BiO layer. This can 
lead to accidental clustering of the dopant oxygen atoms 
which may or may not happen in real materials depend- 
ing on the details of the chemistry. Since very strong 
clustering happens rarely and locally, it does not affect 
the essential physics we describe. 



A. Inhomogeneous electronic structures 

The bare ionic potential calculated from Eq. (Q) is 
shown in Fig. 6a for a 32 x 32 lattice at x = 0.12. The self- 
consistently determined screened potential of Eq. (|lj) is 
plotted in Fig. 6b. The large fluctuations in the bare 
ionic potential are dramatically screened as can be seen 
by comparing Fig. 6a and Fig. 6b. The spatial variations 
of the self-consistently determined LDC in Eq. (Q) and 
Arf in Eq. ( p3| ) are shown as 2D maps in Figs. 7a and 
7b. ft is clearly seen from Fig. 7a that the variations in 
the hole density closely track the screened potential in 
Fig. 6b. The local ionic dopant configurations are shown 
as black dots. The reverse color-coding for the A^-map 
makes it easy to observe the (anti-)correlation between x% 
and Ad(i). The average inter- hole distance at this doping 
is d x ~ 2.86a, which is smaller than the size of the screen- 
ing cloud 2£ x ~ 3.72a of a single isolated ion discussed in 
the previous section. The fact that the number of dopant 
ions is the same as the number of doped holes then im- 
plies that the screening clouds associated with individual 
ions overlap and interfere such that nonlinear screening 
dominates and leads to significant inhomogeneity in the 
electronic structure throughout the entire system. This 
confirms the conjecture we made in Ref. || that the inho- 
mogeneity results from the lack of conventional metallic 
screening of the dopant ionic potential in doped Mott 
insulators. 
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FIG. 7. The spatial variations of the local hole concen- 
tration Xi (a) and the d-wave order parameter A<j(i) (b) 
shown as 2D maps for x — 0.12 on a 32 x 32 lattice. The 
location of the negatively charged ions projected on the 
2D plane is shown as black dots in (a), together with the 
contours of the doped hole density. Reversed color-coding 
is used in (b) for easy visualization of the anti-correlation 
between Xi and Ad(i). 

We present, in Fig. 8, the auto- and cross-correlation 
function analyses for the variations in LDC and the d- 
wave OP defined in Eqs. ([E|[25]J2|) . The 2D hole-hole 
correlation is plotted in Fig. 8a, which shows a rapidly 
decaying central peak and a weak interference pattern at 
larger distances due the presence of many-ionic impuri- 
ties and the d-wave gap structure. Taking the azimuthal 
average of Fig. 8a results in Fig. 8b. The decay length 
extracted from a Gaussian fit is £ x — 1.8a, very close 
to that of the isolated ion case discussed in the previous 
section. This indicates that the main effects controlling 
the average size of the "patches" come from the single 
ion screening cloud. The weak oscillatory structures at 
larger separations is a result of the interference pattern 
in Fig. 8a. The distance to the first weak secondary peak 
can be interpreted as the averaged distance between re- 
gions or patches where the Xi and A<j(i) vary significantly 



from their averaged values. 

The azimuthally averaged auto-correlation for the d- 
wave order parameter variations is plotted in Fig. 8c 
which shows the same structure with a decay length 
^A d — 2.2a. We remark that at this average hole concen- 
tration, which is slightly below optimal doping, the oscil- 
latory structures are rather weak and are just beginning 
to emerge. Interestingly, the azimuthally averaged cor- 
relation functions obtained from the experimental STM 
data on optimally doped BSCCO, shown in Fig. 2 (e-f) 
of Ref. M, have the same structure, including the weak 
oscillations in the tails. We expect the secondary peak to 
develop more prominently in more underdoped systems 
where percolative-like patches are likely to become more 
pronounced with increasing average inter-hole distance. 



correlations indicates that the correlation between the 
spatial variations in A^ and Xi is very strong. 
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FIG. 8. The auto- and cross-correlation analysis of the 
spatial variations in the doping concentration Sxi and the 
d-wave order parameter SAd(i)- (a) 2D auto-correlation 
of Sxi. (b) Azimuthally averaged auto-correlation func- 
tion of Sxi. (c) Azimuthally averaged auto-correlation 
function of 5Ad(i)- (d) The absolute value of the az- 
imuthally averaged cross-correlation function of 8xt and 
8Ad(i)- The solid lines in (b), (c) and (d) are Gaussian 
fits to the correlation functions with the corresponding 
decay lengths £ CT shown in legends. 

Comparison of the ajj-map in Fig. 7a and the Ad-map 
in Fig. 7b shows that the order parameter gap is big- 
ger where the doping rate is lower and vice versa, in 
agreement with the general trend of the doping depen- 
dence of Ad in a clean system. This (anti-)correlated 
behavior is confirmed by the cross-correlation analysis 
of Sxi and 5 Ad in Fig. 8d, where the absolute value of 
the azimuthally averaged correlation function is plotted 
versus separation. The fact that the cross-correlation 
length, Cx-A d — 2.1a, is very close to that of the auto- 
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FIG. 9. (a) A scatter-plot of the d-wave order parame- 
ter Ad(i) and the extended s-wave order parameter A s (i) 
versus local doping concentration Xi. The solid line cor- 
responds to Ad as a function of x obtained in the absence 
of the ionic dopant potential, (b) A histogram showing 
the statistical distribution of Ad(i) in this single 32 x 32 
system. The solid line is a Gaussian fit. 

We proposed, in Ref. |9[, the concept of local doping 
concentration to emphasize the local nature of the physics 
associated with screening and pairing in a doped Mott in- 
sulator. To further substantiate this idea, we produce a 
scatter-plot of the local d-wave order parameter Ad(i) 
versus the local doping xi in Fig. 9a. For a range of dop- 
ing near the average doping value of x — 0.12, the data 
points scatter around the Ad v.s. x curve (solid line) ob- 
tained in the homogeneous case without the ionic dopant 
potential. In the regions where the local doping is low, 
the scatter in Ad is small and the correlation between 
Ad(i) and xi is extremely strong. On the other hand, 
the Ad(i) values show significant scatter in the locally 
overdoped regimes and deviate appreciably from the val- 
ues in a clean system at the corresponding doping level. 
Also shown in Fig. 9a is the scatter plot of the much 
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smaller, extended s-wave order parameter A s (i) versus 
Xi. It scatters around a zero mean and shows a slightly 
larger variance with decreasing LDC. The statistical dis- 
tribution of Arf in this single 32 x 32 system is plotted as 
histograms in Figs. 9b. It can be fitted reasonably well 
by a Gaussian distribution. 



B. Local tunneling density of states as measured by 
STM spectroscopy 

To enable further comparisons with the STM experi- 
ments, we next calculate the local tunneling density of 
states Ni(u;) according to Eq. (pl|). The S functions in 
Eq. ( pl| ) are replaced by the derivatives of the correspond- 
ing Fermi distribution function |p2| at an inverse tem- 
perature of f3J = 30. Parallel to our presentation of the 
experimental data [Bj , we extract the local pairing energy 
gap (full gap), denoted by Ar(i) to differentiate from the 
d-wave OP Ad(i) discussed above, from the peak-to-peak 
distance in the electron LDOS spectrum Ni(cu) at every 
site i. Simultaneously, we obtain the integrated LDOS 
(ILDOS) at every site i from 



U = 



dujNi(uj), 



(27) 



where the integration limit is typically taken to be — luq — 
— J. Note that Ni(u>) at io < corresponds to the LDOS 
of the occupied states below the chemical potential. It 
should be compared to the STM tunneling spectrum at 
negative sample bias where electrons tunnel out of the 
occupied states in the sample. The energy w is related 
to the bias V according to uj = eV . 

The spatial variations of At(i) and 7, are shown as 
the gap and the ILDOS maps in Figs. 10a and 10b re- 
spectively. As used in Ref. [pi, the reverse color coding 
in the gap map (brighter color for a smaller tunneling 
gap) makes it easy to observe its (anti-)correlation with 
the ILDOS map, namely, bigger gap regions correspond 
to smaller ILDOS and vice versa, in agreement with the 
STM data. In Fig. 11a, we construct a scatter plot of 
the superconducting gap versus the ILDOS consisting of 
a total of 32 x 32 = 1024 data points, which suggests an 
approximate linear relationship between the two quan- 
tities, in remarkable agreement with the STM data we 
presented in Fig. 4 of Ref. 0. 

It is easy to notice the great similarity between the IL- 
DOS map in Fig. 10b and the 2D hole density map shown 
in Fig. 7a. This supports within our model the conjec- 
ture made in Ref. M that the ILDOS as measured by 
STM should be approximately proportional to the LDC. 
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FIG. 10. The 2D tunneling gap map (a) and the IL- 
DOS map (b). The reverse color coding in the gap map 
(brighter color for a smaller gap) makes it easy to observe 
its (anti-)correlation with the ILDOS map. The tunneling 
spectra along a line-cut marked by the dashed line and at 
the five locations marked by x in (a) are shown in Fig. 13. 

To examine this relationship in more detail, we show 
a scatter-plot of the ILDOS 7j, obtained from Eq. fl27|), 
versus the local doping level Xi in Fig. lib. Remark- 
ably, the ILDOS is indeed linearly proportional to the 
LDC when x\ is small, i.e. when the local region is un- 
derdoped. This is consistent with the Mott-Hubbard pic- 
ture in which doping a Mott insulator introduces spectral 
weight near the Fermi energy necessary to transform the 
insulator into a superconductor. As observed in recent 
ARPES experiments, [t^ , [l3l]33; 1 the integrated quasipar- 
ticle weight scales with the average doping concentration 
in the underdoped regime. The idea of LDC that we 
advocate extends this Mott-Hubbard picture, commonly 
used to interpret macroscopic spectroscopy properties, 
to microscopic scales. It remains valid in describing the 
local spectroscopies of the short-coherence length super- 
conducting state of the doped Mott-insulator. However, 
for larger hole concentration, i.e. further away from the 
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Mott insulator, the ILDOS in Fig. lib tends to tap- 
per off and saturate. It shows that the two quantities 
are more intricately connected than the overall propor- 
tionality stemming from the quasiparticle wave-function 
renormalization which is apparent in Eq. (pl|). This re- 
sult turns out to be important for the analysis of the low 
energy inhomogeneity structures in the following section. 
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FIG. 11. (a) A scatter-plot of the superconducting gap 
Ar(i) in the tunneling spectrum versus the ILDOS. (b) 
A scatter-plot of the ILDOS versus the LDC x%. 

The statistical properties of the tunneling gap At and 
the ILDOS Ii are in general quite similar to those of the 
d-wave pairing order parameter Ad and the local hole 
concentration Xi. The histogram of At in the 32 x 32 
system is plotted in Fig. 12a, which shows an approxi- 
mate Gaussian distribution. The azimuthally averaged 
auto-correlation of the spatial gap variation SAx(i) is 
shown in Fig. 12b, and that of the cross-correlation of 
the gap variation and the ILDOS variation, 



C,-A r (r,) = ^5>AT««,: + ,) 



(28) 



is shown in Fig. 12c in the absolute value. The extracted 
gap correlation length £a t — 2.9a is somewhat larger 
than that of d-wave order parameter ^A d • Although the 
values of the decay lengths depend on the value of the 
setback distance d s of the dopant oxygen ions, we find 
for d s = 1.5a, and a = 3.8A, £a t — HA which is close 
to the corresponding value of 14A obtained by analyzing 
the experimental data. |9j The similar cross-correlation 



length is consistent with the observation that these two 
quantities are strongly (anti-)correlated. 
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FIG. 12. Statistical properties of the tunneling gap At. 

(a) A histogram of At in the 32 x 32 system at x — 0.12, 
showing approximate Gaussian distribution (solid line). 

(b) The azimuthally averaged auto-correlation function of 
the tunneling gap variations, (c) The absolute value of the 
azimuthally averaged cross-correlation function between 
the variations in At and the ILDOS. The solid lines are 
Gaussian fits. 

Next we turn to the details of the spatial variations 
of the local tunneling spectrum and make comparisons 
to experiments. Fig. 13a shows a three-dimensional ren- 
dering of the tunneling spectra along a horizontal line-cut 
marked in Fig. 10(a), exemplifying the detailed variations 
of the LDOS, the pairing energy gap, and the correlation 
between them. 
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FIG. 13. Spatial variations of the LDOS spectrum as 
a function of bias, (a) A three-dimensional rendering 
of the tunneling spectra along the line-cut indicated in 
Fig. 10(a). (b) The same data from a bird's-eye view. 
The black dotted lines trace the locations of the coher- 
ence peak, (c) Five characteristic spectra taken at the 
positions marked in Fig. 10(a). The local doping concen- 
trations are: Xi = 0.05, X2 = 0.07, X3 = 0.10, X4 — 0.14, 
x 5 = 0.19. 



The projection of Fig. 13a onto the 2D plane is shown 
in Fig. 13b. It is clear that regions with higher ILDOS 
also exhibit smaller pairing gaps and higher coherence 
peaks. The black dotted lines trace the positions of 
the coherence peaks. Combined with the gap-map in 
Fig. 10a, one can see that the gap varies somewhat less 
rapidly near the center of a "patch" than at its "edges" . 
We expect this feature to become more prevalent for 
smaller average doping concentrations, i.e. more under- 
doped samples. In Fig. 13c, we show five characteristic 
tunneling spectra calculated at the five positions marked 
in Fig. 10a. The LDC at these points are given in the 
captions. Notice that the spectral lines plotted versus 
the bias voltage show remarkable resemblance to the ex- 
perimental data shown in Fig. 3 of Ref. ||. The most 
striking feature is that the spatial inhomogeneity exists 
all the way down to low energies as can be clearly seen 
from the the systematic scatter in the zero-bias tunnel- 
ing conductance. [J34| The broad range of agreements be- 
tween the results presented in Figs. 11, 12 and 13, and the 
experimental findings including such details add to our 
confidence that the present theory of the inhomogeneous 
d-wave superconducting state with spinon pairing and lo- 
cal holon condensation captures the essential physics of 
the superconducting state in BSCCO. 



V. ZERO-BIAS TUNNELING DENSITY OF 

STATES 

One of the most striking results of both the STM mea- 
surements (see Fig. 3 of Ref. H) and the present theory 
(see Fig. 13) is that the spatially inhomogeneous elec- 
tronic structure not only appears at the energy scale of 
the superconducting gap, but persists all the way down 
to low energies, including the zero-bias tunneling conduc- 
tance. On the experimental side, the low energy inhomo- 
geneity can only be observed clearly when the tunnel- 
ing matrix element effect along the c-axis is removed by 
a normalization procedure introduced in Ref. J9J. Let's 
take a brief detour and study this matrix element effect. 
The tunneling current at site r^ is given by 



I E (nV) = / duf(u)[l - /(« + eV)] x 

Ni(u)p(w + eV)\M(Ti,z)\ a , 



(29) 



where p(ui) is the density of states of the STM tip which 
is usually taken to be a constant, and Af(r^,z) is the 
tunneling matrix clement. Our knowledge of M is scarce 
in complicated materials, but in general one should be 
able to factor out its dependence on z which corresponds 
to the tip distance to the tunneling surface, 



M(Ti,z) = const, x m(r i )e" Q,5z(r * ) . 



(30) 



Here Szfc) is the spatial variation of the tip-to-surface 
distance in the constant-current scanning mode, and a 
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is a constant determined by the work-function. The de- 
tailed behavior of mfa) is presently unknown and can at 
best be assumed to have no or very weak spatial depen- 
dence. However, the spatial variation in the tip distance 
8z(vi) is known in the form of the surface topography 
from experiments. The normalization procedure we used 
in Ref. || is precisely to remove this part of the matrix 
element effect such that the resulting data represent tun- 
neling measurements carried out with the tip-to-surface 
distance effectively held constant. While it has been com- 
mon practice to present STM data without such a nor- 
malization, in which case the low energy inhomogene- 
ity of the tunneling conductance is not clearly revealed, 
| fiO| , pl| we believe that the constant tip distance normal- 
ization is physical and should be carried out in order to 
obtain the true electronic contribution to the tunneling 
spectra. Our theoretical calculations in the inhomoge- 
neous d-wave superconducting state based on local spinon 
pairing and local holon condensation indeed support the 
experimental observation H that the low energy STM 
tunneling spectra is spatially inhomogeneous. We next 
provide a careful analysis of the theoretical data at zero- 
bias and make a few predictions that can be tested and 
used to guide further experiments. 

Fig. 14a shows a 2D map of the zero-bias LDOS at 
every point on our 32 x 32 lattice. The spatial variations 
are obvious, but what is striking is that this is the same 
pattern of inhomogeneity that we saw at the energy scale 
of the superconducting gap as shown in the gap map in 
Fig. 10a, in complete agreement with experiments. |3q | 
Moreover, comparing Fig. 14a to the local hole density 
map in Fig. 7a shows that they are directly correlated. 
This strongly suggests that the same kind of inhomogene- 
ity pattern should be observable at all energy scales in the 
STM tunneling spectroscopy. To further illustrate the 
correlation between them, we construct a scatter-plot of 
the zero-bias LDOS iV^O) versus the LDC Xi in Fig. 14b. 
It suggests a strong correlation via a linear relationship 
between the two quantities in the regions of low doping. 
The scatter increases gradually with increasing LDC and 
the linear relationship tappers off and saturates. This is 
quite similar to the behavior of the ILDOS and points to 
the more intricate Xi-dependence of the LDOS in Eq. J2l| ) 
beyond the overall linear proportionality. 

Since the local doping concentration is not directly ac- 
cessible experimentally, in Fig. 14c we present the scat- 
ter plot of the zero-bias LDOS versus the ILDOS. In- 
terestingly, due to the weaker variation of the ILDOS in 
regions of higher hole concentration (see Fig. 14b), the 
scatter of the zero-bias LDOS exhibits a slight upward 
curvature with increasing ILDOS. Such a scatter plot can 
be directly constructed from the STM experimental data 
for comparison, [B6[ which serves as another test for the 
present theory. Two remarks are in order. 

(i) It is important to emphasize the difference between 
the electron tunneling DOS that spectroscopy measure- 
ments probe and the thermodynamic DOS or the com- 
pressibility that shows up in thermodynamic measure- 



ments such as the specific heat. While the former de- 
pends on both the quasiparticle wave-function renormal- 
ization Z and the self-energy corrections, the latter is 
insensitive to Z. (21j Although the inhomogeneity in the 
tunneling spectra at low energies emerges from the spa- 
tial variations in the LDC Xi through both the wave- 
function renormalization and the spinon tunneling spec- 
tra as described by Eq. ( pij) , its implications on the ther- 
modynamic DOS and the transport properties require a 
different analysis. We stress that the difference between 
these two density of states even at our unrestricted mean- 
field level is more than the local doping concentration Xi, 
the prefactor in the electron tunneling LDOS, due to the 
presence of the spinon wave- functions in Eq. (Ell) . 
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FIG. 14. Spatial inhomogeneity at zero bias, (a) A 
2D map of the zero bias LDOS, iVi(0), showing the same 
pattern of inhomogeneity as the tunneling gap map and 
the ILDOS map in Fig. 10. (b) A scatter-plot of 7V,(0) 
versus the LDC x t . (c) A scatter-plot of iV;(0) versus the 
ILDOS. 
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(ii) Although the spinon tunneling DOS is not directly 
accessible by STM measurements, since only the phys- 
ical electrons can tunnel in and out of the sample, it 



can be, nevertheless, readily extracted from Eq. (21). 
Dividing out the prefactor x%, we show in Fig. 15a the 
spinon tunneling spectra, N[{ui) = Ni(uj)/Xi, at the five 
locations marked in Fig. 10a. Comparing to the corre- 
sponding spectral lines for electron LDOS at the same 
locations shown in Fig. 13c, it is clear that the degree of 
inhomogeneity remains large at the energy scale of the 
superconducting gap. This is not surprising because the 
gap inhomogeneity results from local spinon pairing in 
our picture. However, the degree of inhomogeneity ap- 
pears to have been somewhat reduced at low energies 
near zero-bias, suggesting that the spinons near the gap 
nodes, unlike the (electron) nodal quasiparticles, expe- 
rience less inhomogeneity. |55| The scatter plot of the 
spinon zero-bias LDOS versus local hole concentration is 
shown Fig. 15b. Although the relative magnitude of the 
variations is reduced by about 30% from the electron case 
in Fig. 14a, the spatial variation of the zero-bias conduc- 
tance for spinons is not only clearly visible, but shows a 
very well defined correlation with the local hole concen- 
tration that is remarkably similar to that of the d-wave 
order parameter Ad shown in Fig. 9a. Therefore, we con- 
clude that the low energy spinons near the d-wave gap 
nodes will experience the same type of inhomogeneity as 
the spinons at the gap edge. 
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FIG. 15. The tunneling LDOS for spinons. (a) The 
spinon tunneling spectra at the five positions marked in 
Fig. 10a, in comparison to those of the electron tunneling 
spectra shown in Fig. 13c where the LDC Xi, i — 1, . . . , 5 
were given, (b) A scatter-plot of the spinon zero-bias 
LDOS versus the LDC Xi, showing remarkable similarity 
to that of Ad(i) versus x% presented in Fig. 9a. 



VI. STM TOPOGRAPHY AND THE 

TUNNELING SPECTRA AT CONSTANT 

CURRENT 



current, which corresponds to the original STM data be- 
fore normalization. In doing so, we will further illustrate 
the physics behind the normalization procedure used in 
Rcf. JEJ, which is essentially a mapping from a constant 
current topography to the electron LDOS measured at 
constant tunneling barrier width. 

The STM measurements are usually carried out at a 
constant tunneling current Jo which is equivalent to the 
integrated LDOS from the sample bias voltage — Vq to 
the Fermi level. Thus we can write, 



Ia = 



dl , 

dV—(r,V,z(r)), 
-v aV 



(31) 



where Jy (r, V, z{r)) is the LDOS or the tunneling differ- 
ential conductance. It is a function of the 2D coordinates 
r on the tunneling surface, the sample bias V, and the tip 
to surface distance z(r) (topography) that must vary at 
every r in order to keep the current Iq a constant. Note 
that J^7 in Eq. ( pl| ) corresponds precisely to the original 

From our discussion 



STM data before normalization 

following Eq. (E 

physical electronic tunneling spectra N(r, eV) obtained 

at constant z(r) according to, 



it is clear that 4rr is related to the 

dV 



N(r,eV) = ^(r,V,z(r)y 



*« 



(32) 



Integratin g E q. ( |32| ) over V from — Vq to gives, after 
using Eq.©, 



7(r) = I e a 



,-(r) 



(33) 



where I(r) is the electronic ILDOS given by Eq. (121 
in terms of the normalized or the calculated tunneling 
spectra. Eq. (pT), together with Eq. (|32"|), completely 
describes the mapping between the constant current to- 
pography and the local electron tunneling spectra at con- 
stant tunneling barrier width. 

From our calculated ILDOS map and using Eq. (|33|), 
we obtain the constant current STM topographic image 
z(r) on our 32 x 32 lattice using a = 2.3A" 1 which corre- 
sponds to a work-function of about 4eV. The topography 
is presented in Fig. 16a, which exhibits identical struc- 
tures as the electron ILDOS map shown in Fig. 10b. We 
next covert the entire local tunneling spectra N(r, V) at 
constant tip-sample distance to the tunneling differential 
conductance at constant current using Eq. ( J32| ) and the 
topography map. The procedure turns out to be equiva- 
lent to normalizing the LDOS spectra by the correspond- 
ing ILDOS. 



We next derive from our numerical data, a theoretical 
STM topographic image at constant tunneling current 
for the 32 x 32 system we study. We then use the re- 
sult to reconstruct the local tunneling spectra at constant 
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FIG. 16. Calculated STM spectroscopy of the 32 x 32 
lattice in the constant current mode, (a) The topographic 
image of our system, showing the same inhomogeneity 
structures as the electron ILDOS map plotted in Fig. 10b. 

(b) The characteristic tunneling differential conductance 
spectra at the same five positions marked in Fig. 10a. The 
integrated area under the spectral lines on the negative 
bias side is the same, i.e. the constant tunneling current. 
Note that the inhomogeneous distribution of tunneling 
differential conductance near zero bias is hardly visible. 

(c) A scatter-plot of the zero-bias conductance versus the 
LDC. 

In Fig. 16b, we show five differential conductance 
curves at the same five positions marked in Fig. 10a. 
Notice that, typical of the raw STM data, the spectral 
lines on the negative bias side cross approximately at 
one value of the voltage as a result of the constraint that 
the integrated area under each curve, i.e. the tunneling 
current, must be the same. The spatial inhomogeneity 



on the energy scale of the superconducting gap remains, 
albeit that the peak height does not show the same sys- 
tematics — smaller gap with higher coherence peak — as 
in the experimental data with or without normalization, 
pa This discrepancy may be due to the simplification of 
the mean-field theory which ignores the effects of fluctu- 
ations, and perhaps also, to a certain extent, due to the 
complication of the matrix element effects in the STM 
tunneling data. The spatial variations of the tunneling 
differential conductance spectra at low energies are, how- 
ever, much less visible. The scatter-plot of the zero-bias 
conductance versus local doping concentration is shown 
in Fig. 16c. Comparing to the same plot for the electron 
zero-bias LDOS shown in Fig. 14b, we find that the rela- 
tive magnitude of the variations in the zero-bias conduc- 
tance is suppressed by about 65% and it shows almost no 
systematic dependence on the LDC. These findings are in 
good qualitative agreement with the experimental obser- 
vation that the original STM data before normalization 
show much less spatial inhomogeneity at low energies. 



VII. SUMMARY AND DISCUSSIONS 

We have investigated in this paper, through micro- 
scopic calculations based on a generalized t-J model, the 
effects of the off-plane ionic potential associated with 
off-stoichiometry doping in the d-wave superconducting 
state of a doped Mott insulator. We find that nonlin- 
ear screening dominates the response of the doped Mott 
insulator to the ionic potential. One of the main char- 
acteristics of nonlinear screening is the emergence of a 
percolative-type of inhomogeneity in the electronic struc- 
ture. In 2D electron systems formed in modulation doped 
semiconductors, nonlinear screening of the ionic dopant 
potential at low electron densities can tear the 2D elec- 
trons into an inhomogeneous mixture of metallic and di- 
electric regions. |14| What we have shown here is a strik- 
ing analogy of the physics in a doped Mott insulator: 
nonlinear screening of the dopant ionic potential leads to 
an inhomogeneous d-wave superconducting state wherein 
the LDC and the d-wave superconducting gap exhibit sig- 
nificant spatially correlated variations on the scale of a 
nanometer comparable to the short coherence length. In 
this context, the concept of a local doping concentration 
augmented by the generalization of the Mott-Hubbard 
picture to local spectroscopies that we advocate || is very 
useful. We have shown that, within our self-consistent, 
spatially unrestricted mean-field approach, local spinon 
pairing and local holon condensation capture the essen- 
tial physics of the inhomogeneous superconducting state 
and provide remarkably consistent descriptions of the ex- 
perimental data. H 

We have shown that there is one length scale controlled 
by the ionic setback distance d s that characterizes the in- 
homogeneity. The decay lengths of the auto-correlations 
of the spatial variations in the LDC and the d-wave pair- 



16 



ing order parameter are found to be close to d s , suggest- 
ing that the superconducting pair-size is determined by 
d s . In this paper, we have chosen d s — 1.5a based on the 
physical distance of ~ 5A between the CuC>2 plane and 
the BiO layer where the dopant oxygen ions reside and 
a ~ 3.8A for the Cu-Cu atomic spacing. However, it is 
important to emphasize the effects of anisotropy in the 
background dielectric screening constant which enters the 
ionic potential. It was pointed out to us by Kivelson that, 
for the high-T c cuprates, the background dielectric con- 
stant, determined by the electronic interband polariza- 
tion and phonon contributions, is highly anisotropic, i.e., 
it has a different value along the c-direction, e±, than in 
the afr-plane e||. If this anisotropy is taken into account, 
the ionic potential given in Eq. (^) must be modified by 
replacing QVi^ V d * and d s — > d* with d* s = ^eJeTd s . 
Thus the effective distance between the BiO layer and 
the Cu02 plane can be significantly larger than 5A. Tak- 
ing e± ~ 8.0 and en ~ 1.5, it follows that the effective 
setback distance d* ~ 3.5a ~ 13A. Our results then im- 
ply that the pair-size, determined by the decay length 
in the spatial variation of the d-wave order parameter, 
would be on the order of one to two nanometers, in good 
agreement with our experimental findings. [|9| 

The physics discussed here can be continued to more 
underdoped cases where the percolative structures result- 
ing from nonlinear screening are expected to be more pro- 
nounced. Already in the case of 12% average doping stud- 
ied here, the tunneling spectra along particular line-cuts 
in Fig. 10a show patches over which the superconducting 
gap is nearly flat in the center and changes quickly near 
the edges, indicative of a percolative electronic structure 
which is often referred to as microscopic phase separa- 
tion. [plnolO| As we have shown that the size of the 
screening cloud around an isolated ion is determined by 
its effective setback distance d*, it is likely that, with in- 
creasing doping, a percolation transition/crossover takes 
place near optimal doping where the average inter-hole 
distance becomes comparable to an effective setback dis- 
tance d*. 

We have not discussed the very underdoped physics 
in detail because in that case it will become important 
to include the effects associated with magnetism that is 
not important and has thus been left out in the average 
doping range considered in this paper. We have tested 
that if magnetism is included, locally antifcrromagnetic 
ordered insulating regions emerge and percolate at suf- 
ficiently low averaged doping. This is left for further 
studies. Nevertheless, the present theory implies that in 
underdoped BSCCO, the local tunneling spectra should 
reveal a percolative structure of superconducting patches 
embedded in a background where larger tunneling gaps 
to single-particle excitations arise. 

A frequently asked question is how does the high 
transition temperature coexist with the inhomogeneity. 
While we do not have a complete answer, we would like 
to discuss a few relevant issues related to this question. 



One of the most obvious reasons for high-T c supercon- 
ductivity to survive the inhomogeneity is its short coher- 
ence length. In our theory, the latter, determined by the 
effective setback distance d* , is also the scale of the in- 
homogeneity and will therefore remain shorter than the 
mean-free path which is usually much larger than the cor- 
relation length of the scattering centers for a reasonable 
scattering strength. Moreover, as we have pointed out 
in Ref. ||, the scattering by the out-of-plane ionic po- 
tential involves predominantly small momentum trans- 
fers limited by the small scattering angle on the order of 
1/dgkp, where kp is the Fermi momentum. This type of 
scattering is much less effective at reducing T c than at 
increasing the single-particle scattering rate. 38 
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FIG. 17. The temperature dependence of the averaged 
d-wave order parameter (open squares) on a 16 x 16 lat- 
tice at an average doping x = 0.15. Also shown in open 
circles is the temperature dependence of the d-wave order 
parameter in the uniform d-wave state without the ionic 
potential. 

We have preliminary results on the temperature depen- 
dence of the distribution of the d-wave order parameter 
Ad- In Fig. 17, we plot the averaged Ad versus tempera- 
ture on a 16 x 16 lattice at an average doping x = 0.15. 
Also shown is the temperature dependence of A^ in the 
same system without the ionic potential. Comparing the 
two curves shows no sign of degradation of T c defined by 
the averaged order parameter. Nevertheless, a distribu- 
tion of Ad does imply a distribution of T c and a finite 
width for the transition. We leave this topic for future 
study. 

There is no evidence at present that would suggest the 
inhomogeneity as being necessary for observing high-T c 
superconductivity. However, the existence of such inho- 
mogeneity should at least help protect the d-wave super- 
conducting state from other possible competing instabil- 
ities. Similarly, from a theoretical perspective, the in- 
homogeneous d-wave superconducting state with spinon 
pairing and local holon condensation may do much bet- 
ter at suppressing low-lying fluctuations beyond the self- 
consistent mean-field solutions than its uniform counter- 
part. 

We next briefly discuss what happens if the dopant 
ions are more ordered in the BiO layer. Experimentally, 
thermal annealing is believed to be able to achieve a cer- 
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tain degree of dopant ordering. To this end, we studied a 
case where the dopant ions are periodically ordered and 
the average doping x = 0.14. We found that the LDC 
and the d-wave order parameter show spatial variations 
that are periodic with the same periodicity as that of the 
dopant ions. The widths of the distribution functions are 
significantly reduced when compared to the disordered 
case. We therefore expect that, although the spatial in- 
homogeneity would remain after thermal annealing, the 
spatial variations would be more periodic in structure 
with a noticeably smaller magnitude. 

Based on these results, we expect that the microscopic 
variations in the electronic properties of YBCO, partic- 
ularly in its ortho-phase, 12] should be more periodic in 
space because the oxygen dopants can be ordered in the 
copper-oxygen chains. However, there exists an impor- 
tant difference in the role of the dopant ions. In BSCCO, 
the BiO plane is essentially insulating which allows us to 
treat the dopants as solely providing the dielectrically 
screened ionic potential. The situation can be very dif- 
ferent in YBCO because the chains themselves are con- 
ducting and have low energy dynamics of their own. The 
electrons in the chains will scatter off the dopant ions 
giving rise to charge density oscillations along the chains 
that may complicate the way by which carriers are doped 
into the Cu02 plane. The present theory would predict 
the appearance of inhomogeneous electronic structures 
in the Cu02 planes of YBCO that are most pronounced 
directly under/above the doped chains. This type of spa- 
tial inhomogeneity on YBCO may have been recently ob- 
served. Hi 

We have investigated here one of the most direct conse- 
quences of off-stoichiometry doping of a Mott insulator, 
namely the ionic potential. There are potentially others. 
For example, the presence of dopants can affect the trans- 
fer integrals in their vicinity, causing spatial variations in 
the parameters t and J, and perhaps more importantly 
in t±_ along the c-axis. It is therefore plausible that cer- 
tain aspects of the inhomogeneous superconducting state 
could also arise from inhomogeneous pairing interactions 
pO| . [Ilf or impurities and defects in the superconducting 
plane W2 43 and many other kinds studied previously 
in the context of disordered superconductors. These are 
interesting and important issues that need to be further 
investigated. 
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